/**
 * This class describes lattice that we use for our Ising model
 */
package simulation;

import java.util.ArrayList;
import java.util.Random;

import common.Neighbors;

/**
 * @author vladimir
 * 
 */
public class Lattice_discrete {
	int nSites = -1;

	public Interval_discrete[][] lattice2;
	public Interval_discrete[] lattice1;
	// int[][] lattice1;

	int[][][][] neighbors2;
	int[][] neighbors1;

	public int coordNumber = 0;
	int dimension = -1;
	double J = 0;
	double B = -1;
	int L = -1;
	double T = -1;
	double beta = -1;
	double LD = -1;
	// int acceptance = 0;

	double epsilon = 0.01;// beta is discretized with this step
	int N = -1;// number of layers

	Random rand = new Random();

	// TupleXY[] sublattice1;
	// TupleXY[] sublattice2;
	// TupleXY[] sublattice3;

	public Lattice_discrete(int[][] neighbors, double J, int L, double B,
			double T) {
		this.coordNumber = neighbors.length;// number of the neighbors
		this.dimension = neighbors[0].length;// dimensionality of the lattice
		this.nSites = (int) Math.pow(L, this.dimension);// total number of
		// sites
		// in spatial direction
		this.L = L;// size of the lattice
		this.J = J;// coulpling constant in spatial direction
		this.T = T;// temperature
		this.beta = 1 / T;
		this.B = B;

		if (B != 0) {
			this.N = (int) (long) Math.max(Math.round(1 / (T * epsilon)), 2);
		} else {
			this.N = 1;
		}
		this.epsilon = 1 / (T * N);

		if (B != 0) {
			this.LD = -0.5 * Math.log(Math.tanh(this.B * this.epsilon));
		} else {
			this.LD = 0;
		}
		System.out.println("prob : " + (1 - Math.exp(-2 * this.LD)));
		System.out.println("LD : " + this.LD);
		System.out.println("N : " + this.N);
		System.out.println("epsilon : " + this.epsilon);

		if (this.dimension == 2) {
			Neighbors NB = new Neighbors(neighbors, this.L, this.L);
			this.neighbors2 = NB.neighbors2;
		} else if (this.dimension == 1) {
			Neighbors NB = new Neighbors(neighbors, this.L);
			this.neighbors1 = NB.neighbors1;
		}

		if (dimension == 2) {
			this.lattice2 = new Interval_discrete[this.L][this.L];
			// create lattice array
			for (int x = 0; x < this.L; x++) {
				for (int y = 0; y < this.L; y++) {
					// fill lattice
					this.lattice2[x][y] = new Interval_discrete(this.B, this.T,
							this.N, this.epsilon, this.LD);
				}
			}
		} else if (dimension == 1) {
			this.lattice1 = new Interval_discrete[this.L];
			for (int x = 0; x < this.L; x++) {
				this.lattice1[x] = new Interval_discrete(this.B, this.T,
						this.N, this.epsilon, this.LD);
			}
		}
	}

	/*
	 * Calculates Energy of the Lattice
	 */
	public double Energy() {
		double result = 0;

		if (this.dimension == 2) {
			for (int x = 0; x < this.L; x++) {
				for (int y = 0; y < this.L; y++) {
					for (int z = 0; z < this.N; z++) {
						result -= this.lattice2[x][y].sites[z]
								* sumNeighbors(x, y, z);
					}
				}
			}
		} else if (this.dimension == 1) {
			for (int x = 0; x < this.L; x++) {
				for (int z = 0; z < this.N; z++) {
					result -= this.lattice1[x].sites[z] * sumNeighbors(x, z);
				}
			}
		}
		result = result * this.J / 2;

		if (B != 0) {
			double temp = 0;
			if (this.dimension == 2) {
				for (int x = 0; x < this.L; x++) {
					for (int y = 0; y < this.L; y++) {
						temp += this.lattice2[x][y].Energy_B();
					}
				}
			} else if (this.dimension == 1) {
				for (int x = 0; x < this.L; x++) {
					temp += this.lattice1[x].Energy_B();
				}
			}

			result += temp;
		}
		return result / this.nSites / this.N;
	}

	/*
	 * Calculates Magnetisation of the lattice
	 */
	public double M() {
		double result = 0;
		if (this.dimension == 2) {
			for (int x = 0; x < this.L; x++) {
				for (int y = 0; y < this.L; y++) {
					result += this.lattice2[x][y].value();
				}
			}
		} else if (this.dimension == 1) {
			for (int x = 0; x < this.L; x++) {
				result += this.lattice1[x].value();
			}
		}
		result = result / this.nSites / this.N;
		if (Math.abs(result) > 1) {
			System.err.println("Magnetisation should be less than 1");
			System.exit(0);
		}
		return result;
	}

	// /*
	// * // * Calculates Magnetisation squared of the lattice //
	// */
	// public double M2() {
	// double Magnetisation = M();
	// return Magnetisation * Magnetisation;
	// }

	// // calculates sum of neighbors within range 2d case
	public double sumNeighbors(int siteX, int siteY, ArrayList<Integer> cluster) {
		double result = 0;
		for (int tn = 0; tn < this.coordNumber; tn++) {
			int neighborX = this.neighbors2[siteX][siteY][tn][0];
			int neighborY = this.neighbors2[siteX][siteY][tn][1];
			result += this.lattice2[neighborX][neighborY].value(cluster);
		}
		return result;
	}

	// calculates sum of neighbors within range 1d case
	public int sumNeighbors(int siteX, ArrayList<Integer> cluster) {
		int result = 0;
		for (int tn = 0; tn < this.coordNumber; tn++) {
			int neighborX = this.neighbors1[siteX][tn];
			result += this.lattice1[neighborX].value(cluster);
		}

		return result;
	}

	// calculates sum of neighbors within range 2d case
	public double sumNeighbors(int siteX, int siteY, int z) {
		double result = 0;
		for (int tn = 0; tn < this.coordNumber; tn++) {
			int neighborX = this.neighbors2[siteX][siteY][tn][0];
			int neighborY = this.neighbors2[siteX][siteY][tn][1];
			result += this.lattice2[neighborX][neighborY].sites[z];
		}
		return result;
	}

	// calculates sum of neighbors within range 1d case
	public int sumNeighbors(int siteX, int z) {
		int result = 0;
		for (int tn = 0; tn < this.coordNumber; tn++) {
			int neighborX = this.neighbors1[siteX][tn];
			result += this.lattice1[neighborX].sites[z];
		}
		return result;
	}

	// updates cluster
	public void update() {
		if (this.dimension == 2) {
			for (int x = 0; x < this.L; x++) {
				for (int y = 0; y < this.L; y++) {
					for (int z = 0; z < this.N; z++) {
						int tx = rand.nextInt(this.L);
						int ty = rand.nextInt(this.L);
						ArrayList<Integer> newCluster = this.lattice2[tx][ty]
								.newCluster();

						double dS = DeltaS(tx, ty, newCluster);

						if (rand.nextFloat() < Math.exp(-dS)) {
							this.lattice2[tx][ty].set(newCluster);
						}
					}
				}
			}
		} else if (this.dimension == 1) {
			for (int x = 0; x < this.L; x++) {
				for (int z = 0; z < this.N; z++) {
					// for (int z = 0; z < this.N; z++) {
					int tx = rand.nextInt(this.L);
					// int tz = rand.nextInt(this.N);
					ArrayList<Integer> newCluster = this.lattice1[tx]
							.newCluster();
					// double dS = DeltaS(tx, tz);
					double dS = DeltaS(tx, newCluster);
					// System.out.println("dS : " + dS);
					if (rand.nextFloat() < 1 / (1 + Math.exp(dS))) {
						this.lattice1[tx].set(newCluster);
						// this.acceptance += newCluster.size();
						// this.lattice1[tx][tz] *= -1;// .sites[tz] *= -1;
						// }
					}
				}
			}
		}
	}

	// calculates change in Energy between old and new configuration
	public double DeltaS(int x, int y, ArrayList<Integer> newCluster) {
		int tempValue = this.lattice2[x][y].value(newCluster);
		if (Math.abs(tempValue) != newCluster.size()) {
			System.err.println("in DeltaS, cluster is not uniform");
			System.exit(0);
		}
		// return this.epsilon * 2 * this.J * tempValue
		// * sumNeighbors(x, y, newCluster);
		return this.epsilon * 2 * this.J * tempValue
				* sumNeighbors(x, y, newCluster) + 2
				* this.lattice2[x][y].DeltaE_B(newCluster);
	}

	public double DeltaS(int x, ArrayList<Integer> newCluster) {
		int tempValue = this.lattice1[x].value(newCluster);
		if (Math.abs(tempValue) != newCluster.size()) {
			System.err.println("in DeltaS, cluster is not uniform");
			System.exit(0);
		}
		// return this.epsilon * 2 * this.J * tempValue
		// * sumNeighbors(x, newCluster);
		return this.epsilon * 2 * this.J * tempValue
				* sumNeighbors(x, newCluster) + 2
				* this.lattice1[x].DeltaE_B(newCluster);
	}
}
